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Abstract 

This is the second of two papers which investigate cold dilute neutron matter on the lattice 
using pionless effective field theory. In the unitary limit, where the effective range is zero and 
scattering length is infinite, simple scaling relations relate thermodynamic functions at different 
temperatures. When the second virial coefficient is properly tuned, we find that the lattice results 
obey these scaling relations. We compute the energy per particle, pressure, spin susceptibility, 
dineutron correlation function, and an upper bound for the superfluid critical temperature. 
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I. INTRODUCTION 



This is the second of two papers which investigate cold dilute neutron matter on the 
lattice using pionless effective field theory. We refer to the first paper as [I]. In this paper 
we study the scaling behavior of the results in the unitary limit and compute a number 
of physical observables, the energy and pressure, the spin susceptibility, and the dineutron 
correlation function. 

In the unitary limit the scattering length is much larger than the interparticle spac- 
ing whereas the range of the interaction is much smaller. In this limit the only length 
scales in the problem are (n^ (/i))~ 1//3 and the thermal wavelength Ay. Here, nf(p) = 



(37r 2 ) -1 (2mjv/x) 3 / 2 is the density of a free Fermi gas and At = a/2 71 "/ (m N T), where is 
the neutron mass, /i is the chemical potential, and T is the temperature. All dimensionful 
quantities can be expressed as suitable powers of either (n/( / u)) -1 / 3 or Ay times a function 
of the dimensionless quantity p/T. For example, we can write the pressure as Q 

2 

P{T,n) = -nn f {fi)g{x) (1) 

where x = fi/T. Using standard thermodynamic identities one can show that 

3P 



p(p) = n f (n) 



(2) 



where p is the density and e is the energy density. In this work we first show how these 
relations arise in the lattice regularized theory. We then present numerical results for the 
equation of state, the spin susceptibility, and the dineutron correlation function, and show 
that the lattice data are consistent with universality. 



II. SCALING IN THE UNITARY LIMIT 



In this section we derive some scaling relations which relate observables at different tem- 
peratures in the unitary limit. Our derivation is equivalent to the treatment in For 
the analysis it is simplest to first take the temporal lattice spacing a t — > 0, and then take 
the spatial lattice spacing a — > 0. When we take the temporal lattice spacing a t — > 0, we 
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end up with a Hamiltonian lattice formulation, 

n s ,i 

~2^Y1 Yl \4{^s)ai{n + i s ) + a\(n s )ai{n s -i s 

ft,i Z s =l,2,3 

+ a\(H s )a^(n s )a[(n s )ai(n s ). (3) 

Here, a>i(n s ) is an annihilation operator for a neutron with spin index i at the spatial lattice 
site n s and C is the coupling constant. 

In the unitary limit, which corresponds with a — > and a sca tt — > oo, it is not difficult to 
show that 

C=~5-, (4) 

m N 

where r\ is a constant, 

L 3 

= lim — ~ 3.957, (5) 

°° integer f2- 

27rfci 27r/C2 27rfc3 

= 6 — 2 cos 2 cos 2 cos . (6) 

k L L L w 

The derivation is as follows. At zero temperature, the value of C can be set by the condition 

that the two-particle scattering pole occurs at the energy prescribed by Liischer's formula 

for energy levels in a periodic box of length L In the actual lattice simulations we 

describe later we will determine C in a slightly different way, but the two procedures agree 

in the limit that the lattice spacing goes to zero. If we place the two-particle scattering 

pole at energy E po \ e then the condition on C is 

~C = L-^L L3 2 _ E ■ ( 7 ) 

k integer P m " k 

Liischer's formula gives 

£pole = ^^[l + 0(^)], (8) 

We keep a sca tt finite for the moment. In the limit L — > oo we split the sum in (J7J) into the 
term k = and the remaining terms k ^ 0, 

1 11 1 v ^ rriN 

k=^0 integer 



m N , m N ^ 1 



47ra srat t L 3 ^ fL 



fc^O integer k 
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Therefore 



C 



m N 



1 • 



47ra scat t L% Z^fe^O integer 

k 



In the unitary limit where a sca tt — > oo, we have 

1 L 3 



C 



m N Ek^O integer ft m JV ' 



which corresponds with the three-dimensional attractive Hubbard model 

\U\ 

^ = 2ri ~ 7.914. 
t 

The grand canonical partition function for our system is 



(10) 



(11) 



with 



(12) 



Z G = Trexp \—(3(H — fiN)] 



(13) 



or 



Zq = Tr exp 



En, (^-^)al(n) ai (n) 



En,i Ej=i,a,3 a l(n)ai{n + I) + a\(n) ai (n - 1) 

Eh a \ (™W ( H ) a \ ( H ) a i («) 



We observe that is a function of only two parameters, and /3/i. In terms of physical 
units, these parameters are 

3 i i / , . \ - 



(14) 



2m w 2mP) lys TP h y s a 2 4vr ( At ' ' 



(15) 



phys 



yphys 



In z, 



where z is the fugacity, 



z — e 



Bp, 



(16) 



(17) 



We can take the two independent parameters to be Ay hys a _1 and z. 

Consider any operator F(a{, aj) built from the lattice annihilation and creation operators. 
We assume that 



continuum 



lim 

a->0 



a D ■ ( F(as, a 



(18) 



has a nonzero continuum limit for some power D. We know that 



F(a i ,al 



Tr 



F(oi,a\)exp[-P(H-iiN)) 



Tr [exp {-(3{H -fiN)]] 
4 



(19) 



is a function of only Ay hys a 1 and z. So in order that the factor of a D drops out of (jl8|). we 
need in the continuum limit 

(F( ai ,a]))^(\^ s a- l ) D f(z), (20) 
where f(z) is some function of the fugacity only. Hence 

(^)cont_=(A? yS ) D /W (21) 

and so 

(^T 1 ^) (^continuum ( 22 ) 

is a function of z only. We extend this to the case when the operator has an explicit 
dependence on the displacement, x 9 ^. In that case 

<^(- PhyS )>co nt „ - S k • (*K 4 ^a- 1 ))} , (23) 



where x^^a 1 is the displacement in lattice units. Therefore 

(A?-) ~° <^(^ hys )> conti _ = / (z, ^(Arr 1 ) • (24) 

As an example we show that the particle density times three powers of the thermal 
wavelength is a function of only the fugacity. The particle density in the continuum limit is 



(p) continuum = [a ■ {a\ (n s )a T (n s ) + aj(n s )a^(n s ) \J , (25) 

where the expectation value can be measured at any spatial lattice site n s . Therefore 
D = — 3 and 

^continuum (26) 

is a function of only the fugacity. Similarly the energy per particle times inverse temperature, 

£§; (27) 

pressure times inverse temperature and three powers of the thermal wavelength, 

\v hys Yf3P; (28) 



and Fermi energy times inverse temperature, 
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are all functions of fugacity only. 

Using the unitary limit scaling relations we can derive a simple relation between the 
energy density and pressure The energy density is 

E 1 d 



V 



Vdp 



In Z G + fip 



1 

V 



d n d 

dp ~ pajl 



\nZ G . 



We note that 



3 1 
— In Zn 
V 



,phys 



which is only a function of fugacity. Since 

d n d 
~dp ~ pdfl 

we have the result 



d /j, d 
dp ~ fid/l 



/3P, 



J" = 0, 



E 
V 



1 

V 



d yU d 

dp ~ pdjl 



\nZr 



(30) 



(31) 



(32) 



~ te^Viiu z< 



V 
3 1 

2pV 



InZ 



V 
3 



d n d 
~d~p ~ pdjl 



\phye 



G 



P. 



(33) 



III. RESULTS 



In the following we present lattice simulation results for cold dilute neutron matter in 
the unitary limit. We use a spatial lattice spacing of a = (50 MeV) -1 and temporal lattice 
spacing of a t = (24 MeV) -1 . The temporal lattice spacing is sufficiently small that the 
results are close to the a t — > limit. We use the hybrid Monte Carlo algorithm j?J to 
generate Hubbard- Stratonovich field configurations as described in |8|. We use diagonal 
preconditioning before each conjugate gradient solve as described in [I]. 

The finite volume error was tested by going to larger volumes, and the final lattice sizes 
were chosen so that the finite volume error was less than one percent. The lattice dimensions 
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we used for the various temperatures are shown in Table 1. 

Table 1: Lattice volumes for various temperatures 



T (MeV) 


L 


L t 


4 


4 


6 


3 


5 


8 


2 


5 


12 


1.5 


6 


16 


1 


6 


24 



We adjust the coefficient C of the four-fermion interaction as suggested in [I]. For each 
temperature and lattice volume, we compute the density for both free lattice fermions and for 
lattice regularized two-particle bubbles. We use the free fermion result in order to determine 
the first virial coefficient b\(T) (see Table 1 in [I]) and the bubble sum to compute the second 
virial coefficient b 2 (T), 

(34) 



bubble 



T3-&1CO [z + 2.b 2 (T)z 2 ]. 

Arp 



We then adjust C so that 



b 2 (T) = 3-2-3. 



(35) 



This constraint is also used to compute the derivative of C with respect to the temporal 



lattice spacing. The results for C(T) and f=- are shown in Table 2. 



dat 



Table 2: C and ^ on the lattice 



T (MeV) 


C (lO" 4 MeV" 2 ) 


(10~ 5 MeV- 2 ) 


4 


-0.971 


-0.76 


3 


-0.948 


-1.40 


2 


-0.958 


-2.12 


1.5 


-0.987 


-2.47 


1 


-1.043 


-2.68 


0.667 


-1.098 


-2.65 


0.5 


-1.128 


-2.52 



For each temperature we have probed densities up to a quarter-filled lattice. With a spatial 
lattice spacing of (50 MeV) -1 , the quarter-filled lattice corresponds with a density of 0.0081 
fm -3 . Beyond this one might find significant lattice artifacts. 
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FIG. 1: Density times A^ versus fugacity for various temperatures. 
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FIG. 2: Density times versus fugacity for small fugacity. We also include comparisons with 
the virial expansion at first and second orders. 



A. Density versus fugacity 

In Fig. ^we plot the density times A^ versus fugacity for temperatures T = 4, 3, 2, 1.5, 
and 1 MeV. The data from the five different temperatures appear to fall on a single curve. 
This suggests that p\\ depends only on z, as predicted by unitary limit scaling. In Fig. 
121 we magnify the plot of pA^ at low fugacity, showing both bubble chain results and full 
simulation results. We show the first order and second order virial results with 
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FIG. 3: Energy per particle times f3 versus fugacity. We show results for T = 4, 3, 2, 1.5, and 1 
MeV. 



b 2 (T) = 3-2 _ § ^ 0.530. (36) 

Since we have tuned the interaction coefficient to produce the correct second order virial 
coefficient, it is not surprising that the lattice data agrees with the virial expansion at low 
fugacity. 



B. Energy per particle versus fugacity 

In Fig. El we plot the energy per particle times (3 versus fugacity for temperatures T = 4, 
3, 2, 1.5, and 1 MeV. The energy per particle times (3 appears to depends only on fugacity, 
as predicted by scaling in the unitary limit. The small deviation for different temperatures 
appears to be due mainly to an overall shift in the height of the curves. In the continuum 
limit at z = 0, we expect the equipartition result 

¥-§■ 

The slight deviations from | at z = can be attributed to lattice cutoff effects in the free 
particle kinetic energy. 

In Fig. 0] we magnify the same plot for small fugacity and include both bubble chain 
calculation results and full simulation results. We also show the first and second order virial 
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FIG. 4: Energy per particle times (3 versus fugacity at small fugacity. We compare with first and 
second order virial expansion results. 



results in the continuum. At first order we have 

PE 3 



A - 2' < 38) 



and at second order, 

PE_ (§-M* + ^(f-21n*)** 

~t ~ : T 3 „ 2 + ln z - ^ 

We see that apart from small shifts in the overall height, the lattice results agree with the 
continuum virial results. 



C. Energy density and pressure 

In Fig. El we show the energy density times (3\% versus fugacity for temperatures T = 4, 
3, 2, 1.5, and 1 MeV. Scaling in the unitary limit requires that the energy density times 
(3\t is only a function of fugacity, and this appears to be the case. In Fig. |S]we also plot 
the pressure times |/3A^, which according to ()33|) should equal the energy density times 
(3\t- We have computed the pressure by numerical integration of the density as a function 
of chemical potential, 

P= ^lnZ G = 1 T A(ji')dfj! = pQi'W- (40) 

* ' J — CO J —CO 
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FIG. 5: Energy density times (3\p versus fugacity for various temperatures. We compare with the 
pressure times |/3A^. 



We see that the lattice results appear to confirm the unitary limit relation, 



D. Reduced energy versus reduced temperature 

In units where Boltzmann's constant equals 1, the degeneracy temperature Tp is the same 
as the Fermi energy, 

Tf . Ef .W!. (42) 
2m N V ; 

In Fig. El we plot the energy per particle divided by \Ep versus the temperature divided by 
Tp for temperatures T = 4, 3, 2, 1.5, and 1 MeV. As expected from unitary limit scaling 
all points appear to lie on a single curve. In Fig. [7|we show a magnified plot of the energy 
per particle divided by \Ep versus the temperature divided by T F . The data points at the 
lowest values of T /Tp appear to lie on straight line with intercept 

= 0.07. (43) 

One expects the curve to reach T/T F = with zero slope. The actual intercept at T/T F = 0, 
which corresponds with the parameter £ in the zero temperature relation, 

A 5 2m 
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FIG. 6: The energy per particle divided by \Ep versus temperature divided by Tp for temperatures 
T = 4, 3, 2, 1.5, and 1 MeV. 
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FIG. 7: A magnified plot at low temperatures of the energy per particle divided by \Ep versus 



temperature divided by Tp. 



should likely be somewhere between 0.07 and 0.42. Since there is no noticeable nonanalytic 
behavior in the data shown, this suggests that the superfluid critical temperature Tq should 
satisfy 

^<0.14. (45) 

T F 
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FIG. 8: Lattice results for the axial response as a function of fugacity for temperatures T = 4, 3, 
2, 1.5, and 1 MeV. 



E. Axial response and spin susceptibility 



The axial response at zero momentum is defined as 



91 



Sa(0) = -L m _ Nim _ ATJ) , (46) 

where iV^jj is the number operator for spin up(down) neutrons. is the total number 
operator so that 

(N) = (JV T + iV ; ) = A. (47) 

S'a(O) is a dimensionless quantity with a finite continuum limit and therefore should be a 
function of fugacity only. It is normalized to equal 1 at zero density. Sa(0) is proportional 
to the Pauli spin susceptibility j^, |j| , 

Xp = E ( [al(n s )[a) ijaj (n s )] • [aKfifMkMK)] ) = (48) 

rt s ,n' s 

In Fig. [HI we show the axial response as a function of fugacity for temperatures T = 4, 3, 
2, 1.5, and 1 MeV. The lattice results for the axial response appear to lie on one curve, as 
predicted by the unitary limit scaling relations. 

The decrease in the axial response indicates that the transfer of spin from one spatial 
region to another is being suppressed. This is due to the formation of spin zero pairs. A 
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FIG. 9: Logarithm of the dineutron correlation function versus lattice distance measured along a 
lattice axis. We show results for T = 1 MeV. 



superfluid transition leads to nonanalytic behavior in the susceptibility. A sharp crossover 
in the susceptibility can be used to define a pseudogap phase. We do not observe any of 
these phenomena. The susceptibility becomes quite small for the most degenerate systems 
studied, but the dependence on temperature is smooth. 



F. Dineutron correlation function 

We define the equal-time dineutron correlation function as 

G#(n 8 ) = (^ai{n s )a^{n s )a\{Q)a\{^)j . (49) 

In Fig. El we show the logarithm of the dineutron correlation function as a function of 
lattice distance measured along a lattice axis. We show results for T = 1 MeV and various 
values of T/Tp. We have staggered the plots for better viewing, and the five highest points 
represent data measured at zero lattice distance. For comparison we also show the free 
dineutron correlation function as well as the result of the bubble chain approximation. We 
observe that the interacting correlation function is larger than the non-interacting one for 
all temperatures. This reflects the attractive interaction in the spin zero channel. At the 
higher temperatures the correlation function is quantitatively explained by the bubble chain 
approximation. At T/Tp = 0.14 the bubble chain correlator starts to level off, indicative of 
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long range order, but the full correlation function does not. This implies that Tc/Tp < 0.14. 



IV. SUMMARY 

We have studied neutron matter in the unitary limit on the lattice. In [I] we concentrated 
on the low density/high temperature equation of state and compared the results to the 
virial expansion. We found significant finite lattice spacing errors and suggested that the 
scaling behavior can be improved by tuning the second virial coefficient rather than the zero 
temperature scattering length. 

In part II we showed that tuning the second virial coefficient improves universal scaling 
for the entire range of densities and temperatures studied. We found, in particular, that 
the energy per particle in units of Tp only depends on fi/T and that the energy density is 
3/2 times the pressure. In the temperature range studied the energy per particle in units of 
Tp is essentially linear in T /Tp. This means that at present we can only provide bounds on 
the universal parameter £. We find 0.07 < £ < 0.42. This result is marginally consistent 
with the Green function Monte Carlo result £ = 0.44 |lo| . 

We also studied the spin susceptibility and the dineutron correlation function. We find 
that the spin susceptibility is strongly suppressed in the most degenerate systems, but no 
clear phase transition is observed. We also find that the dineutron correlation function is 
strongly enhanced over the free correlation function, but no long range order is observed. 
We conclude that the critical temperature for the the transition to a superfluid phase is less 



than 0.14Tp. This bound is lower than the result Tq = 0.22(3)Tp 



but consistent with 



the result Tq = 0.035(4)Tp obtained in Q|. Green's function Monte Carlo calculations are 
restricted to T = and can only determine the gap, not the critical temperature. Carlson, 
et. al., [ipj find A = 0.9 • \E F . If the critical temperature were related to the gap by 
Tc = 0.57A as in BCS theory this would imply Tc = 0.31Tp, but there is no reason to 
expect this relation to hold in the unitary limit. 
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